1 Introduction 



Evolution of massive stars (M >~ 20M) is accompanied by a mass loss, initiated 
by a high luminosity and high radiation pressure. In blue supergiants, situated 
near the main sequence, mass loss rate is moderate M ~ 10~ 6 MQ/yr. and 
is connected with the outflow of layers having small optical depth. This mass 
loss rate is determined by radiation pressure in lines, where spectral absorption 
coefficient may be very high. Theory of the outflow is based hear on Sobolev 
approximation of the lines treatment, ( CAK theory ), and velocity gradient in 
the flow increase essentially the pushing effect of the radiation pressure (Castor, 
Abbot, Klein , 1975). First attempt to construct a theory of such mass outflow 
have been done by , and for recent development of CAK theory see (lO) . 

Evolved massive stars may lose mass with much higher rate then the blue 
supergiants. Formation of single Wolf-Rayet stars probably took place as a result 
of such intense mass loss (Bisnovatyi-Kogan, Nadyozhin, 1972). Outflowing 
atmospheres of luminous massive stars are described by equations of radiation 
hydrodynamics with radiative heat conductivity. The flux goes through a critical 
point where the velocity is equal to the local isothermal speed of the sound. 
The main difference between the outflowing atmospheres in blue and evolved 
supergiants is a value of the optical depth ( average Rosseland ) in the critical 
point - r cr . In blue supergiants we have r cr < 1, and acceleration is essential 
only in parts of the radiation spectrum around the spectral lines, but in evolved 
yellow supergiants there is r cr >> 1 and the whole radiation flux takes part in 
the flow acceleration. 

Theory of the mass outflow from the very luminous evolved stars was devel- 
oped by Bisnovatyi-Kogan, Zeldovich, (1968), and self - consistent method for 
evolutional calculations with a mass loss was proposed by (Bisnovatyi - Kogan, 
Nadyozhin, 1969, Bisnovatyi - Kogan, Nadyozhin, 1972). The method was based 
on finding the solution of the outflowing envelope which continuously passes the 
isothermal speed singular point and smoothly enters the second singular point 
at infinity where v — > t>oo, T — > 0. The shortcoming of the consideration of 
[pi, where evolution with self - consistent mass loss was first calculated, was 
ignorance of the fact, that the optical depth in the flow is decreasing and far 
from the star r — ► there is no influence of matter on the radiation. Formally 
equations used in this paper described only optically thick outflow and could 
not be extended to infinity. 

The aim of this paper is to correct this shortcoming and to give a unified 
description of the flow which is started at large optical depth and goes to infinity 
at r = 0. The approximate system of equations based on Eddington approxi- 
mation for the radiation is derived and solution is found which is continuous in 
both singular points. 

The recipe based on the approximate treatment of the outer boundary con- 
dition was proposed for the problem of the mass loss by |la| , [n6| . Similarly to 
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the case of a static atmosphere it was suggested that at the stellar photosphere 
where T — T e ff = (i/47rr 2 cr) 1 / 4 there is a fixed value of the optical depth 
T p h = 2/3, or later r P h = 8/3 in the papers of (Kato, 1985, Kato, Iben, 1992). 
The basic approximation made in these papers was based on the approximate 

oc 

choice of the value of r p h , which instead of its expression r = J npdr was taken 

r 

as Tph ~ Tph = {npr) p h- It looks out here, that we avoid all problems connected 
with the singular point at infinity, because the region of solution is restricted 
by r < r p h- In fact this approximation creates another problems. In the case 
of incomplete ionization and rapid change of the opacity as a function of T and 
p this approximation could introduce rather big errors into the solution. The 
method considered is free from such shortcomings also. In the solution pre- 
sented below we have T p h = 4 and f p h — 3.75. Recently another approximate 
approach for solution of the outflowing atmosphere was proposed by |L3| Ex- 
panding layers of the star are divided into subsonicaly extending photosphere 
where the stationary momentum equation is adopted, and the wind where the 
accelerating up to velocity of the outflow is approximated by the prescribed 
profile with 3 parameters which are determined from matching and additional 
conditions. 



2 Thermodynamic relations in outflowing enve- 
lope at arbitrary optical depth. 

In deep layers, at large optical depth r>l we have a usual equation of state 
of a mixture of ideal gas with a black - body radiation. 

P(p,T) = ^+pKT. (1) 

When t is small it is not correct to speak about equilibrium between radiation 
and matter. Since local thermodynamical equilibrium is still presumed we take 
isotropic radiation component together with a gas pressure in the form 

nT 4 

P™f = — • r + pKT. (2) 

Multiplier r in (2) ensures zero decreasing contribution of a component at 
the infinity. 

An increasing anisotropy of a pressure of radiation coming from the star at 
decreasing t is taken into account in Eddington approximation. Basing on the 
solution of the radiative transfer equation in Eddington approximation for the 
spherically symmetric case with variable Eddington factor (Sobolev, 1967), we 
introduce the following approximate representation of thermodynamic relations 
for radiation 



2 



H __ T , , L(r) 



^- — (1-0 + ^ (3) 

P £V = «r*(i-c-) + ^. (4) 

Term in (3), (4) containing L t h is to be determined from a solution of self- 
consistent set of equations . This term essentially contributes only when r is 
small. This fact allows us to simplify further calculations by writing out 
instead of L t h{r). 

3 Basic equations. 

A system of equations of radiation hydrodynamics describing continuous transi- 
tion between optically thick and optically thin regions for the stationary outflow 
is written as: 

du _ 1 dP g GM(l - L th ) 
dr p dr r z 

, f Lth(r) T AncGM 
where L th = — , L ed = 

Led K 



L = ^(s+?--^ + ^j+L th (r) (6) 

4nr 2 c/dP r E r p-3P r \ 

Lth = —, (7j 

up \ dr r J 



M 

47T 



p = pur 2 (8) 

„rp4 T OO 

T OO 

Ep = aT ^l-e-^ + ^t c +E g (10) 

Pg = PUT (11) 

E g = l-RT (12) 



■ / 



npdr (13) 



Where L - is a constant total energy flux consisting of the radiative energy 
transfer together with the energy of the matter current, u is a rate of the outflow, 
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k is an opacity, assumed to be constant, a is the constant of a radiative energy 
density, 1Z is a gas constant. 

We consider here the flow in the gravitational field of constant mass M, neglect- 
ing self-gravity of the outflowing envelope outside the critical point m en . This 
approximation is very good in a realistic case of m en << M (Zytkow, 1973). 
This system of equations provides a description of a stationary outflowing enve- 
lope accelerated by a radiative force at arbitrary optical depth where continuum 
opacity is prevailed. In optically thick limit r — > oo, when terms in (9) , (10) 
with are negligible, and E r p = 3P r a solution of this system was obtained by 
[0. In case of small r from (9), (10) for the anisotropic radiation flux we get: 
E r p ~ P r , what follows from transfer equation in Eddington approximation. 
(Sobolev, 1967). 

When optical depth becomes small separation of radiation and matter should 
be taken into consideration. It means that only a part of radiation is interacting 
with outflowing gas envelope. For this part of quantums we still assume LTE 
to be valid, what means that we take the temperature in (3) - (4) to be the 
same as in (11) - (12). For the rest part of radiation that is not interacting with 
outflowing matter, another 'temperature'- mean energy of quantums should be 
introduced. This part of radiation transfers momentum to the outflowing matter 
and thus produces only the anisotropic part of the pressure, that is of the order of 
Lth/r 2 . We avoid the problem of treating this another 'temperature' by solving 
the radiative transfer equation in momentum form (7). It is solved together 
with the equation of motion (5), energy conservation (6) and thermodynamical 
relations for the outflowing gaseous matter and radiation (9) - (10). 

Equation of motion (5) is written in the form where the effect of the pres- 
sure gradient is explicitly written only for the gas component and all effects 
of radiation is embedded into the last term. This term is obtained from the 
transfer equation (7) when taking into account relations (9), (10). In case of 
large optical depths this approach insures that acceleration is providing by the 
total (isotropic = radiation + gas ) pressure gradient, and in case of small r by 
gas pressure gradient and pressure of radiation flux. In intermediate region this 
terms will compete in accordance with (3), (4) and (7). About the difference of 
this system of equations from previously used see the Discussion. 

Substituting (7) - (10) into (6) we get 
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(14) 



Differentiating in (14) with taking into account that dr = —npdr will result in 
a first differential equation 
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i 2n -tJT J 5 GM p? 
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t- ^ar 2 cT 4 e~ T . 
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Here we introduce a coefficient of a heat conductivity 

4acT 3 

A 



3ftp 



(15) 



(16) 



fn a limiting case of r — ► oo this equation coincides with a corresponding 
equation from Q , when taking into account that 
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Anr 2 



< caT 4 



Using (f 1) , (8) in (5) will get 



(TIT /i 2 \ rfp rfT 2 2 GM - 

H 7 = _K T + ' J T5 — (1-Ltfc). (17) 

\ p p°r^ / dr dr p z r° r z 



The system (|T^) , ( |T7| ) need to be completed from (13) by relation 



dr 

— = -Kp. 

dr 



(18) 



The equation ( |17| ) has a singular point, where the left hand side of it vanishes 

KT cr p 2 



Per Ptr r cr 



(19) 



This point corresponds to the "isothermal sonic" point where 



Using (U) , (7) gives 



L th = -4?rr 2 



<9P S 
dp ) T 



A^(l-e~ r ) - -acT^e- r 
dr 3 



(20) 
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The second singular point of ([l5|) - ( plj| ) is situated at infinity, where 

T = 0, p ^ — —> 0, u ~ > const — Wqo, when r — > oo. (21) 

This condition is related to the fact that far from the star the density in stellar 
wind is very small. In reality the wind may be treated as stationary only up to 
the limiting radius r; im ~ v^t >> r cr , where t is the characteristic mass- loss 
time of the star. So the formal solution with outer boundary condition (21) is 
very close to the real solution with conditions at ru m . This approximation is a 
common way to consider a well-developed solar wind problem (Parker, 1963). 
Let us introduce nondimcnsional variables 

T(r) = » , p(r) = £> , L th = ^, (22) 

± rr Oct 



^cr 
X = 

r 



After transformations we obtain dimensionless system of equations 



1 A ie - T T 

A 3 + 7 A +2^7 

4 A 5 x 2 p 2 



dT (h , lx A , _ . . T 4 

^ = -{2 T - A * x+ 27 + {1 - e )Al T (24) 



^ T T4 ' 2L-A 3 A 5 X - - ^ 



AA 2 A 5 x 2 ° p A 2 )T3(1~ 

dr p 



(25) 

Where L°° = L^f t . To simplify writing here and further we omit tilde. Dimen- 
sionless coefficients Ai are: 

Al -3^ ' ^"l^I*' (26) 
, GM ^ = ^J!- (27) 



r cr lZT cr ' 16ac7r r cr T^ r 
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Physical sense of At parameters have been revealed in jjj . Additional parameter 
is 

A 5 = —L- , (28) 

T cr Kp cr 

the reciprocal 1/A 5 is of the order of optical depth in the critical point. 
Condition of transition of the solution through the critical point confines the 
number of dimensionless parameters. Equating to zero expression in the figure 
brackets of (23) in the critical point with account of (24) will get 



A 4 = (4A 3 A 5 (1 - e" T -) + 8A 5 (-1 + e~ T ^) 

+4AlA 2 A 5 (l - 2e- T <" + (e^-) 2 ) (29) 
+A 2 {12A 5 + A 3 (-4A 5 + 8A 2 5 L°°)) 
+AiA 2 (16 A 5 (l - e- T °-) 
+A 3 (4A 5 (-l + e- T -)+A 2 5 (8L 00 
-8e-^L°°)))+e- T "-)/4(A 5 + A X A 5 {1 - e~ T -)). 

All dimensional parameters of the flux: T , p , r could be expressed as a func- 
tion of dimensionless parameters Ai and a dimensional combination of physical 
constants IpJ 



4an 2/5 (GM) 7 / 5 1 
3c ' ^8/5 {AlA 2 A 3 Y^A 3 x' 



55 ^z^^" ,^ 

' V 4a \kGM J K 1 2 S> A l 



In numerical calculations coefficient ^5 should be expressed as a function of A\ 
- A 3 and a nondimensional combination of physical constants and thus it is not 
an independent parameter 

/o\l/5 / 4 1 /5 7? 4/5 

4> = ( 7 J s7s 4/. ■ ( 3 °) 

4 Numerical solution 

In order to satisfy boundary conditions far from the star we need to integrate 
(23) - (25) from the critical point outward to the infinity. We exit the critical 
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point by means of expansion formula. Expanding the solution in critical point 
x = T = p = 1 in powers of (1 — x) we have 

T=l + a(l-x), (31) 



p= 1 + 6(1 -x). (32) 

Similarly 

e- T -e- T 11 + |), 

where y = 1 — x. 

For the a and b coefficients we get 



b = (-12A 2 A 5 - iA 1 A 2 A 5 + AA 2 A 3 A 5 + 4A 4 A 5 

+AA 1 A 2 A b e- T " r - SA2A3AIL 00 (33) 
-e—)/(4A 5 (-l + e-^)), 
a = (- Cl - (c? - 4c 2 co) 1/2 )/(2c 2 ). (34) 

The coefficients q due to their complicated form are adduced in Appendix. 

Numerically integrating we escape the critical point by means of the expan- 
sion formulas |3l]), ( p2| ) . Then integrating outward to the infinity we satisfy 
the boundary conditions . Results of the numerical calculations are shown 
in Fig. 1-2. Curves at this figures correspond to the following dimensional values 
of the dimensionless parameters: A\ = 50, A 2 — 10~ 4 , A3 = 43.88, r cr = 125, 
= 0.6. This set of parameters corresponds to the following values at critical 
point: T cr = 1.4 • 10 i K, r cr = 2.6 ■ 10 13 cm, p cr = 6.6 ■ lO^g/cm 3 . 

The behavior of the solution with rapidly decreasing Mach number at r < r cr 
shows, that it may be matched to the static solution for the core. In reality the 
opacity peak is situated near critical point, opacity inside is decreasing and the 
inside velocity drop is more rapid, then in the case of re = const (Bisnovatyi 
- Kogan and Nadyozhin, 1972). The effective temperature of the photosphere 
may be obtained from: L^/Aur 2 — erT 4 . For the given set of parameters we 
get: x p h — 0.03, r p h — 4, T e ff, = 0.06, and f p h = 3.75. Conditions in critical 
points impose two relations upon the set of nondimensional parameters. As it 
was shown the solution depends upon the following nodimensional parameters: 
A%, A 2 , A3, A4, r cr , Lf^. When numerically treated A4 was expressed as a 
function of the remaining dimensionless parameters when taking into account 
that in sonic point x = p = T = 1 ((2^). Parameter A3 is to be numerically 
determined from the condition at the infinity where T = (pl|). 

When r — » 00 velocity tends to constant and thus p ~ 1/r 2 . Only the 
unique value of T cr allows to obtain the proper behavior of u (or p) at the 
infinity. Varying r cr , we adjust the behavior of u to get u — > Uoo at r — > 00. 
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At r = oo the solution could be represented in an expansion form . In this 
case could be directly determined from (|20|), and the appropriate expansion. 

At the stellar core optical depth may be taken arbitrary large, and it is 
necessary to match only T and p. Above mentioned method of obtaining the 
unique value of r cr allows to avoid the problem of fitting r when integrating to 
the stellar core. 

All other dimensionless parameters of the envelope and parameters of the 
static core could be obtained only when matching the solution for the stellar core 
with the solution for the expanding layers. In this approach all the treatment 
will be fully self - consistent. 

In our treatment we have not used the second expansion and thus could not 
specify uniquely the L^. 

Our aim was to develop a method of a calculating the parameters of the flux 
in arbitrary r. Fully self - consistent treatment which may be applied to the 
real star is under the consideration. 

5 Discussion 

Solution for the spherically-symmetrical stationary outflowing envelope acceler- 
ated by the radiative force in arbitrary optical depth case was obtained in this 
work. 

We have introduced thermodynamical relations for the matter flux with par- 
tially separated radiation. This approach provides satisfactory description of 
the problem for the arbitrary r. In case of r — > oo equation of state yields 
E r p = 3P r and our method reproduces results of Also when r — > 0: 

E r p ~ P r accords the result obtained from the radiative transfer equation in 
Eddington approximation. 

As a result of this treatment we have introduced a system of differential 
equations. This system provides a continuous transition of the solution between 
optically thin and optically thick regions. To satisfy boundary conditions the 
solution should proceed through the critical point where the speed of the flux 
equals the local isothermal sound speed. We have derived analytically approxi- 
mate representation of the solution at the vicinity of the sonic point. Using this 
representation we numerically integrate the system of equation from the critical 
point to satisfy conditions at the infinity. 

Beginning from (Zhytkow, 1973), the obtaining of outer boundary conditions 
was oversimplified, making a wrong impression, that regions with r < 1 , have 
no influence on the mass outflow. Zhytkow imposed boundary conditions at 
the photosphere. That was also the result of that there was no method for 
the self-consistent description of radiation and matter at small t. It should 
be mentioned that most of the papers only inherited the method developed by 
Zhytkow. As soon as it is possible to describe correctly the whole region of the 
flow one should forthwith impose correct boundary conditions. Problem of the 



infinite boundary conditions was elaborated for the theory of Solar wind. 

In papers of Zytkow, it was made an attempt to describe the whole flow. 
The description was based on Paczynsky approximation (Paczynski, 1969), that 
rouphly takes into account the dilution of the radiation flux. For r < 2/3 for 
the part of the equation of motion, that describes pressure, it was taken: 

dp dP a dT up L , dT L 1 „. ._ _i/o o /9 

£ ~ -w* + £^ where * - -i^A - 2 f{T)T ° R ° r ' ' 

where Tq = LjA^acR\ - the effective temperature of the photosphere, and L 
is constant for r < 2/3, / = 1 — (3/2)r. The second term in the relation for 
the temperature gradient originated from the approximation of regions of small 
t ( with Eddington factor 1 ), with regions of t >> 1 ( with factor 1/3 ) and 
is to describe roughly the process of the radiation dilution. The shortcoming 
descended from the description of Zhytkow concerns this last effective surface ( 
with T = T e ff, t = 2/3 ). The 'temperature' of free, not interacting quantums 
corresponds to the prescribed T e ff of the photosphere. It can be not a bad 
approximation, if mostly neglecting the influence of the regions with small r, but 
to describe the whole flow, this approach seems to be too rude. It seems obvious 
that the solution is artificially restricted by taking this prescribed separation 
with the certain value of the effective temperature of the photosphere. 

When r becomes small, separation of the radiation and matter progresses 
in two parts (two 'temperatures', see appropriate discussion), with accordance 
to the transfer equation in momentum form (7). If not treating the transfer 
equation, together with equations of hydrodynamics, the energy conservation 
will be broken for the radiation at the regions of small r. That will only enlarge 
the uncertainty within the extended atmosphere. 

The main difference of our approach with previously published is that we 
use the equation of motion, in which the effects of the gas pressure are sepa- 
rated from the effects of the radiation, together with the transfer equation (7). 
The effect of the gas pressure gradient that is valid at small r as well as in 
deep interior is explicitly written. At r » 1 the main driving force originates 
from the gradient of the pressure of the gas together with the radiation. Far 
from the star acceleration occurs mainly due to the momentum transfer from 
the radiation to matter. In the way the equation (5) is written it is correct for 
arbitrary optical depth, if, of course, to calculate L th from (7). This description 
provides the proper competition of this essentially different effects of the radia- 
tion. Contrary to other authors which used approximations obtained from (7), 
we are resulting from the 'equations of state' (9) - (10). We believe, that this 
approach is almost free from the shortcomings mentioned above. 

We have used equation of state which corresponds to constant R g . Account 
of variable R g , M and k are of principal importance for the real stars as it was 
described by (Bisnovatyi-Kogan, Nadyozin, 1972). This complication does not 
change our method and qualitative results. When one will take into account 
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lines effects, the importance to describe correctly the region of small r will 
forthwith be grater. 

As it was mentioned at the beginning, the problem of taking f instead of 
t needs to be treated more carefully. In this paper we simply assumed k to 
be constant. When considered far from the star velocity of the flux may be 
approximately taken constant, and thus from M = Anpur 2 , get p ~ 1/r 2 . In 
this case r ~ f. Otherwise if we consider a power law of k = p a it is easy to 
obtain r ph = (npr) ph /(\ + 2a), what is (1 + 2a) times smaller than f p h- For 
steeper (exponential) decrease of k with r the difference between r p h and f p h 
may be much larger. 

In real stars opacity may increase more steep. Even if temperature is de- 
creasing smoothly, partial recombination of ions will cause an increase of k ( 
Cox, Tabor 1975 ). Hence it will create regions where flow is accelerated and u 
can not be taken constant and thus law for p will be far from 1/r 2 . In this case 
rather big errors could be introduced when taken f instead of r. 

On a supergiant phase massive stars may have regions of not fully ionized 
H and He what cause an increase of the opacity which leads to the acceleration 
in continuum ( Bisnovatyi - Kogan, Nadyozhin 1969, 1972 ). 

To take into account real physical effects our method should be improved 
by taking non-constant k and variable ionization degree. If improved this way 
theoretical approach presented in this work should be useful for self - consistent 
simulations of the massive stars evolution with mass loss. 

This work was partly supported by RFBR grant 96-02-17231, grant 96-02-16553, 
and CRDF grant RP1-173, Astronomical Program 1.2.6.5 . 

References 

[1] Bisnovatyi-Kogan G.S., 1967, Prikl. Mat. Mech., 31, 762 

[2] Bisnovatyi-Kogan G.S., Nadyozhin D.K. 1972, Ap&SS, 15, 353 

[3] Bisnovatyi-Kogan G.S., Nadyozhin D.K. 1969, Nauch. Inform., 11, 27 

[4] Bisnovatyi-Kogan G.S., Zeldovich Ya.B., 1968, AZh, 45, 241 

[5] Castor J., Abbot D., Klein R., 1975, ApJ, 195, 157 

[6] Cox A.,Tabor J., 1976, ApJS, 31, 271 

[7] Kato M., Iben I., 1992. ApJ, 394, 879 

[8] Kato M., 1985, PAS J, 37, 19 

[9] Lucy L., Solomon P., 1970, ApJ, 159, 879 

[10] Owoki S., Pols J., 1996. ApJ, 462, 894 



11 



[11] Paczynski B, 1969., Acta Astron., 19, 1 

[12] Parker E.N., Interplanetary dynamical processes. Interscience Publishers, 
New York - London, 1963 

[13] Shaerer D., do Koter A., Schmutz W., Maeder A., 1996, A&A, 310, 837 

[14] Sobolev, V.V., Course of theoretical astrophysics. Moscow, Science, 1967 

[15] Zytkow A., 1972., Acta Astron., 22, 103 

[16] Zytkow A., 1973., Acta Astron., 23,121 



12 



6 Appendix 

Further sophisticated calculations were made using Mathematica 2.2 . 
c = (Ai^A§Af96(-l + e- T ^){L°°) 2 + A3(-40A2A 3 i°° + A 2 A 2 A 3 (-160 + 320e- T -- 

160(e- T -) 2 )i°° + A 2 (A 3 176(l - e - r -)i°° + A 2 (-48 + 48e- T -)L 0O ) + 

A 1 (A 2 A 3 (80 + 96A 4 + (-160 - 96A 4 )e- T - + 80( e - T -) 2 )i oc + A|(A§96(1 - er^L 00 + A 3 392(-l - 
2( e -^=r)2 + A 5 (2(e- T ='') 2 + e- T -(-2 + 8A 4 ) + A 2 e- Tc ''(-29 + 8A 3 ) + 2A 3 (l e - T ^ - (e- T ") 2 )+ 
2AlA 2 e- Tcr (l - 2(e~ rCT ) + (e^-) 2 ) + Ai(-2Aie- T "- + 2A 4 (e- TCT ) 2 + A 2 (-7e- T ^ + 7(e" Tc -) 2 + A 3 
{^ Tar ) 2 )))) + -4 2 (-96 + Aj(-60 + 20A 3 ) - 72A A + (192 + 72A A )e~ T ^ - 96(e~ T -) 2 + + 24A 4 - 

(-48 - 24A 4 )e- T - + 24( e - r -) 2 ) + A\Al(-h6 + 168e- T - - 168(e- T -) 2 + 56(e" T -) 3 )+ 
A\{A 2 {2A + 80A 4 + (-72 - 160A 4 )e- T - + (72 + 80At)(e- T -) 2 - 24( e - T -) 3 ) + Aj(-292 + 584e- T - 
A 3 (80 - 160e- T - + 80(e- T -) 2 ))) + A 2 (232 + 20A 4 + A§(24 - 24 e - T -) - 232 e - T -+ 
A 3 (-152 + 152e- T - - IQer^L™)) + A X (-2AA± - 24^1 + (48A 4 + 24A 2 >- T "' - 24A 4 (e- T -) 2 + 
AK-392 + A 3 (196 - 196e- T -) + 392e- T - + A§(-24 + 24e- T -)) + A 2 (160 + 196A 4 + (-320 - 196,4 
160(e- r -) 2 + A 3 (-56 - 48A 4 + (112 + 48A 4 ) e - r - - 56( e - T <-) 2 + L°°(4e- T ^ - 4( e - T ") 2 )))))) /(2(- 
ci = Al(8 + A2(56 - 20,4 3 )- 20A 4 + A 3 (12-12e- r -)-8e- T - + 
A\A 2 {% - 16e- T - + 8(e- T -) 2 )+ 

A X {-\2A± + 12A 4 e- T "- + A 2 (48 - 48e- T <- + A 3 (-12 + 12 e - T -)))) + 
A\ {32A 2 A 3 L°° + A^Asim - IQe-^L 00 ) + 5A 5 e- T -, 
c 2 = A 2 (8 - 8e- T ") . 
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7 Firure captions 



Fig.l Results of the integration of equations (23) - (25) from the critical point 
x = 1 outward to the infinity x — 0. Points x = 1 and x = are the singular 
points of our system. A solution shown for A\ = 50, = 10~ 4 , A3 = 43.89, 
r cr = 125, = 0.6, passes through the critical point (T cr = 1.4 • 10 4 -ftf, 
r cr = 2.6 • 10 13 cm, p cr = 6.6 • 10~ 12 .g/cm 3 ) and satisfies boundary conditions 



Fig. 2 Dimensionless velocity (solid line) and Mach number = u ■ 



(dashed line). Inwardly decreasing Mach number ensures matching the static 
core even in the case of k = const. When r — > 00 velocity tends to constant 
~ 16km/ s and p ~ 1/r 2 . Solution is shown for the same values of nondi- 
mentional parameters as on Fig. 1. 
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